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Abstract:  We  present  a  survey  of  results  from  an  extended  project  focused  on  the 
understanding  of  the  dynamic  behavior  of  elastomers  or  filled  rubbers.  This  en¬ 
tailed  experimental,  modeling,  computational  and  theoretical  efforts.  Of  particular 
emphasis  are  the  nonlinear  and  hysteretic  aspects  of  dynamic  deformations. 
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1  Introduction 

In  this  presentation  we  summarize  joint  research  efforts  between  scientists  at  the  Thomas 
Lord  Research  Center  of  the  Lord  Corporation  and  applied  mathematicians  in  the  Center 
for  Research  in  Scientific  Computation  (CRSC)  at  N.C.  State  University  beginning  in  1994 
and  continuing  through  the  current  efforts  by  the  authors  of  this  report.  The  research  par¬ 
ticipants  from  Lord  were  Lynn  Yanyo,  Mike  Gaitens,  Beth  Munoz,  and  Oon  Hock  Yeoh, 
while  significant  CRSC  contributors  include  H.T.  Banks  (1994-  ),  Yue  Zhang  (1994-1997), 
Nancy  Lybeck  (1995-1997),  Laura  Potter  (1996-1998),  Gabriella  Pinter  (1997-  )  and  Negash 
Medhin  (2000-  ).  As  should  be  apparent  from  our  outline  below  (presented  in  a  somewhat 

1  htbanks@eos.ncsu.edu 

2  Corresponding  author 

3Plenary  Lecture,  The  47th  European  Study  Group  with  Industry  and  the  Mathematics  for  Industry 
Workshop,  August  24-29,  2003,  Graasten,  Denmark. 

4ngmedhin@eos. ncsu.edu 

5gapinter@uwm.edu 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2003 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2003  to  00-00-2003 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Multiscale  Considerations  in  Modeling  of  Nonlinear  Elastomers 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

North  Carolina  State  University, Center  for  Research  in  Scientific 
Computation, Raleigh, NC, 27695-8205 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBIECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

17 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Multiscale  Considerations 


2 


chronological  format),  this  project  is  a  prime  example  of  one  wherein  a  data-driven  applica¬ 
tion  subsequently  leads  to  new  computational,  theoretical  and  conceptual  ideas.  It  results 
in  a  multiscale  framework  that  connects  pseudo-phenomenological  or  phenomenological  ap¬ 
proaches  to  molecular  considerations.  Our  focus  in  these  efforts  has  been  the  understanding 
of  nonlinearity  and  hysteresis  and  the  roles  they  play  in  the  dynamic  deformations  of  filled 
rubber. 

An  important  area  where  nonlinear  models  have  a  potential  impact  is  the  rubber  indus¬ 
try.  Elastomeric  materials  can  be  found  in  diverse  engineering  applications,  e.g.,  in  springs, 
bearings,  shock  absorbing  bushes,  helicopter  rotor  suspensions,  tires,  vibration  suppression 
devices,  bridge  supports.  They  appear  both  as  passive  damping  devices  and  actively  control¬ 
lable  smart  materials.  These  new  applications  motivate  the  need  for  a  better  understanding 
of  the  mechanical  behavior  of  rubber-like  composites  which  is  a  necessary  first  step  in  the 
design  of  both  passive  and  active  material  devices  [1,  2],  The  dynamic  behavior  of  these 
elastomers  is  very  complex.  They  exhibit  significant  nonlinearities  in  both  geometric  and 
material  characteristics.  Typical  nonlinear  behaviors  of  the  stress  and  strain  in  rubber  ma¬ 
terials  under  finite  deformation  include  a  continuous  increase  of  strain  at  decreasing  rates 
upon  loading,  variable  magnitudes  of  strain  subject  to  rates  of  loading  and  different  loading 
and  unloading  paths  due  to  hysteretic  memory  effects.  In  addition,  the  current  state  of  ma¬ 
terial  depends  on  the  strain/strain-rate  history,  the  type  and  amount  of  filler  in  the  material 
and  the  temperature.  The  nonlinear  effects  are  particularly  strong  for  large  strains  and  for 
highly-filled  elastomers.  We  have  attempted  to  accurately  capture  the  nonlinear  dynamic  as 
well  as  hysteretic  effects  which  we  discuss  briefly  next. 

Hysteresis  is  a  widespread  phenomenon  in  science  and  engineering.  It  plays  significant 
roles  in  electromagnetics  (polarization,  conductivity)  in  dielectric  and  conductive  materi¬ 
als,  in  biological  systems  (time  delays  in  disease  pathogenesis,  intracellular  metabolic  path¬ 
ways,  tissue  viscoelasticity),  in  recent  investigations  of  ecological  migrations  (diffusion  with 
integro-differential  terms  representing  local  episodic  behavior)  as  well  as  in  the  deformations 
of  polymeric  materials  such  as  filled  rubbers.  Hysteresis,  often  referred  to  as  “memory”, 
does  not  involve  physical  memory  mechanisms,  but  rather  is  a  manifestation  of  “hidden”  or 
unmodeled  internal  or  local  dynamics.  It  is  an  embodiment  of  the  multiscale  aspects  of  mod¬ 
els  and  phenomena.  That  is,  one  encounters  Preisach  to  Boltzmann  to  internal  dynamics  to 
molecular-based  models  as  one  moves  from  macro  to  micro  to  nano  and  molecular  formula¬ 
tions.  Thus,  approaches  range  from  the  phenomenological  involving  an  input/output  view¬ 
point  to  a  pseudo-phenomenological  involving  internal  dynamics  formulations  of  Boltzmann 
hysteresis  operators  to  physics-based  models  (e.g.,  the  stick-slip  models  of  Doi/Edwards  and 
Johnson/Stacer  that  are  discussed  below)  at  the  molecular  level  for  internal  dynamics. 

2  Models  for  elastomeric  materials 

Traditionally  there  are  two  approaches  to  modeling  rubber  materials.  There  are  physics  or 
molecular-based  theories  [3,  4,  5]  that  try  to  describe  the  microscopic  behavior  of  the  particles 
and  fibers  that  constitute  the  material  and  there  are  phenomenological  models  that  treat 
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the  material  as  a  continuum  [6,  7,  8,  9,  10,  11].  Phenomenological  models  are  based  on  the 
assumption  that  the  material  is  isotropic,  i.e.,  the  long  chain  molecules  are  randomly  oriented 
in  the  unstrained  state.  Loading  causes  orientation  of  the  molecules,  but  it  is  in  the  direction 
of  the  loading,  so  the  assumption  of  isotropy  remains  valid.  Most  such  models  utilize  a  strain 
energy  density  to  describe  the  state  of  the  material.  However,  these  strain-energy  functions 
can  typically  capture  only  the  current  state  of  the  rubber  and  cannot  distinguish  between 
different  loading  and  unloading  histories.  Thus  these  models  cannot  accurately  describe  the 
significant  hysteresis  exhibited  by  filled  rubber  samples. 


2.1  Model  development 


In  developing  a  model  for  the  dynamic  behavior  of  filled  elastomers  we  considered  two  basic 
deformations:  extension  and  shear.  Our  hysteretic  models  were  based  on  the  basic  models 
developed  by  Banks  et.  al.  in  [12,  13].  First  we  considered  the  model  of  a  rubber  rod  under 
uniaxial  tension.  (Our  shear  results  are  summarized  in  Section  2.4.)  The  Timoshenko  theory 
for  longitudinal  vibrations  of  a  rubber  bar  with  a  tip  mass  leads  to 
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(3) 


where  u  is  the  deformation  in  the  x  direction,  p  is  the  mass  density,  f{t )  is  the  applied 
external  force,  Ac  is  the  cross-sectional  area,  M  is  the  tip  mass  and  g  is  the  gravitational 
constant.  This  model  includes  a  Kelvin- Voigt  type  term  as  a  first  approximation  to  damping. 
The  stress-strain  relationship  in  the  basic  model  is 


<r(i)  =  G(e(t),e(t)), 


where  e  is  the  finite  strain  (since  we  are  interested  in  large  deformations)  and  it  is  given  by 


du  1  /du\  2 
dx  2  \dx) 


=  £  + 
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where  e  =  is  the  usual  infinitesimal  strain  of  linear  elasticity.  However,  modeling  the 
nonlinear  behavior  between  the  stress  and  the  finite  strains  e  (which  are  themselves  nonlinear 
functions  of  the  infinitesimal  strains  e)  can  be  equivalently  formulated  in  terms  of  nonlinear 
relationships  between  the  stress  and  the  infinitesimal  strains  e.  Hence,  we  have  equivalently 
assumed  in  obtaining  (1)  a  stress-strain  relationship 


.  .  .du  d2u  .du.  d2u 

a{t)=9’{di’didi)=9’{di)+CD — 


dtdx 


(4) 
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Figure  1:  Elastomer  rod  under  tension  for  dynamic  free  release  experiments. 


instead  of  a  similar  law  involving  the  finite  strains  e. 

The  well-posedness  of  the  initial-boundary  value  problem  associated  with  (l)-(3)  was 
shown  by  Banks,  Gilliam,  Shubov  in  [14],  while  the  convergence  properties  of  the  corre¬ 
sponding  parameter  estimation  problem  were  established  by  Banks,  Pinter  in  [15].  These 
results  are  valid  for  a  broad  class  of  nonlinearities.  This  is  important  since  comparison 
between  experimental  data  and  numerical  calculations  suggested  that  the  neo-Hookean  non¬ 
linearity  (found  in  literature  as  a  first  approximation  to  the  nonlinearity  exhibited  by  rubber 
materials)  is  not  adequate  to  describe  the  behavior  of  filled  elastomers  [13]. 

An  adequate  form  of  ge  and  the  unknown  constants  p  and  Cp  were  determined  using 
parameter  estimation  techniques.  Data  for  the  inverse  problem  were  provided  by  dynamic 
free  release  experiments.  The  elastomer  was  suspended  vertically  with  the  top  end  (x  —  0) 
fixed,  and  a  frame  was  attached  to  the  lower  end  (see  Figure  1).  Varying  amounts  of  extra 
mass  were  attached  to  this  frame,  which  also  served  to  house  an  accelerometer.  Another 
accelerometer  placed  at  the  top  of  the  sample  was  used  to  verify  the  clamped  boundary 
condition  at  the  fixed  end.  For  the  free  release  experiment,  the  rubber  rod  was  lifted  together 
with  the  frame  and  the  tip  mass  so  that  no  compression  or  extension  occurred.  Then  the 
support  was  removed,  allowing  the  mass  to  fall  freely.  This  type  of  experiment  was  repeated 
with  unfilled  and  lightly  filled  carbon  black  samples,  while  a  similar  experiment  was  done 
with  a  highly  filled  sample  with  a  9.29  lb  tip  mass. 

The  force  data  collected  by  a  load  cell  on  top  of  the  sample  were  used  in  estimating  the 
unknown  parameters  q  —  ( ge,p,Cn )  in  (l)-(3)  by  minimizing 

1  N 

j{q)  =  \z*  ~  0;  q)  |2 

1  i= l 

over  q  in  some  admissible  parameter  space  Q.  Here,  z,,  i  =  1 . . .  N,  represent  the  experimental 
observations  at  times  U  of  the  force  at  the  fixed  end,  and  a  is  given  by  (4),  where  u  is  the 
solution  of  (l)-(3)  corresponding  to  the  parameters  q. 
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Computational  results  indicated  that  (i)  a  nonlinear  function  ge  is  necessary  in  the  stress- 
strain  relationship,  (ii)  the  stress-strain  relationship  (4)  is  adequate  to  describe  unfilled 
rubber  samples  (Figure  2)  and  (iii)  to  capture  the  dynamics  of  filled  elastomers,  hysteresis 
must  be  taken  into  account.  Based  on  models  for  hysteretic  damping  in  the  literature  (see 


Figure  2:  (a)  Time  domain  approximation  and  (b)  the  FFT  of  the  solution  vs.  the  data: 
Model  with  a  four-term  piecewise  linear  ge. 

[16]  and  references  therein)  we  assumed  a  Boltzmann-type  stress-strain  law  of  the  form 

<r{t)  =  ge(s{t))  +  CDi(t)  +  J  Y(t- s)-^gv(s(s),e(s))ds,  (5) 

where  Y  is  the  memory  kernel,  and  ge  and  gv  are  nonlinear  functions  accounting  for  the 
elastic  and  viscoelastic  nonlinear  responses  of  the  elastomer,  respectively.  This  stress-strain 
law  implies  that  the  stress  depends  not  only  on  the  current  strain  and  strain  rate  but  also  on 
the  history  of  the  strain  and  the  strain-rate.  It  is  very  important  to  note  that  the  stress-strain 
law  (5)  contains  various  standard  internal  strain  or  internal  variable  formulations  as  special 
cases.  The  anelastic  displacement  field  (ADF)  models  of  Lesieutre  [17,  18]  for  composite 
materials  exhibiting  both  elastic  and  anelastic  displacement  fields  are  formulated  on  the 
assumption  that  the  host  elastic  material  contains  anelastic  materials  with  internal  strains 
£\  which  are  elastic  strain  driven.  That  is,  the  constitutive  laws  have  the  form 

a(t)  =  Eie(t)  -  E2£i{t), 

where  the  internal  strain  is  given  by 

£l(t)  +  ~£i(t)  =  C2£(t),  £i(0)  =  0,  (6) 

T 


£i 


c2e 


~£(s)ds. 


or  equivalently, 
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Several  generalizations  of  this  formulation  exist,  e.g.,  Johnson  et.al.  [19],  suggest  that  the 
internal  strain  is  strain  rate  driven,  i.e. , 

ii(t)  +  -£i(t)  =  c2i{t).  (7) 

T 

Our  Boltzmann-type  law  (5)  (under  appropriate  assumptions  on  the  past  memory  from  —  oo 
to  0)  corresponds  to  an  internal  strain  model  of  the  form 

ii(t)  +  ^£i(t)  =  e(t)),  ei(0)  =  0.  (8) 

This  form  was  chosen  after  we  found  that  neither  (6)  nor  (7)  provided  laws  that  could 
describe  our  data. 

Experimental  observations  of  the  quasi-static  behavior  of  elastomers  indicate  that  these 
materials  possess  different  nonlinear  viscoelastic  responses  in  loading  (i  >  0)  and  unloading 
(i  <  0).  This  led  to  our  choice  of  a  piecewise  continuous  form  for  the  viscoelastic  response 
function  gv 


gv(e(s),i(s)) 


gvi(e(s ))  e{s)  >  0 
gvd(e(s ))  e(s)  <  0. 


(9) 


Here  we  require  gvi  and  gvd  to  be  continuous  (and  generally  nonlinear)  functions.  The 
difference  between  loading  and  unloading  is  more  pronounced  as  the  amount  of  filler  increases 
in  the  material.  We  define  points  f,,  i  —  1, . . . ,  M,  as  the  “turning  points,”  or  the  points  in 
time  for  which  e  =  0.  The  function  gv  need  not  be  continuous  at  the  turning  points,  so  we 
must  interpret  the  derivatives  in  (5), (8)  as  distributional  derivatives.  That  is,  delta  functions 
are  involved  in  differentiating  the  composite  functions  gv(e,i),  or  equivalently,  integration 
by  parts  is  valid.  For  experimental  and  computational  purposes,  we  further  assume  that 
any  motion  in  the  material  prior  to  each  experiment  is  negligible.  That  is,  we  assume  that 
s{t))  ~  0  for  all  t  <  0.  Hence  we  may  approximate  (5)  by 


<7 


IS. 


r(f)  =  ge(e(t))  +  CDe(t)  +  [  Y(t  -  s)-£-gv(e(s),e(s))ds 

Jo  as 

We  integrated  by  parts  in  (10)  and  obtained  the  model  in  variational  form 

,du ,  ,  d2u 


(10) 


„d2u  d 

p-r^-A, 


dt2 


dx 


f  ,du .  .  .  ,du  d2u  .  rt 

S‘(di)+  + 


r  •  ,  .  ,du  .  .  6  u  ,  ..  , 

Lnt-s)9-(As)'mFx{s))is 


+0-V) j,,;.,  +  ^  r(t  —  tk)(— j— (ft))  —  9vd{ —  F{t)  in  r*  (11) 

(12) 


du 

u(t,  0)  =  0,  u(0,a:)  =  v?o,  =  0 


for  an  appropriately  chosen  Hilbert  space  V.  This  presumes,  of  course,  that  we  have  sufficient 
smoothness  so  that  evaluation  of  at  U  makes  sense  and  -§^(gvi( ^(U))),  ~§^{9vd(^(U)))  € 
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V*.  For  this  particular  example,  one  takes  V  =  0,£)  =  {0  £  L2(0,£)\cf)'  £  L2(0,£),  0(0)  = 

0}  and  p  =  pAc  —  M5 ^  F(t)  =  where  Se(x)  is  the  Dirac  function  with  atom  at 

x  —  £.  The  important  question  of  well-posedness  of  the  system  (11)-(12)  is  treated  carefully 
in  [20].  We  have  shown  that  under  certain  assumptions  there  exists  a  unique  weak  solution 
of  the  system  (11)-(12).  This  result  uses  the  Minty-Browder  monotonicity  method  to  treat 
the  nonlinearities  in  the  system  and  is  valid  for  a  broad  range  of  elastic  and  viscoelastic 
response  functions. 


2.2  Experimental  and  computational  results 

We  first  tested  the  hysteretic  model  (11)-(12)  in  the  quasi-static  behavior  of  filled  rubber 
samples  under  uniaxial  tension.  Two  different  types  of  quasi-static  pull  tests  were  conducted. 
The  first  type  (Type  I)  includes  a  sequence  of  loading  and  unloading  the  sample  to  produce 
load-displacement  curves  with  decreasing  maximum  strain  levels.  In  the  second  type  of  ex¬ 
periment  (Type  II)  we  created  a  sequence  of  strain  loops  that  have  decreasing  maximum 
strain  levels  as  before,  but  instead  of  a  fixed  minimum  strain  level  we  used  progressively 
increasing  minimum  strain  levels.  (Full  description  of  the  quasi-static  and  dynamic  tests 
and  results  are  given  in  our  paper  [21].)  We  tried  a  number  of  linear  and  nonlinear  func¬ 
tions  for  ge  and  gv,  including  the  special  cases  of  ge  and  gv  linear  with  gv  =  gvi  =  gvd. 
The  relative  errors  discussed  in  [22]  suggested  that  nonlinear  functions  are  necessary  for 
both  ge  and  gv.  Additional  curve  fitting  studies  outlined  in  [22]  led  to  our  choice  of  cubic 
polynomials  for  ge,gvi  and  gvd.  That  is,  we  chose  parameterized  nonlinearities  of  the  form 
9e{x)  =  E,3=i  Eix\  gvi(x )  =  Y,Uiaix\  9vd(x )  =  E-=i  bix\  where  Ehaubu  i  =  1,2,3  are 
real  constants.  Note  that  these  response  functions  include  no  constant  terms.  Here  we  re¬ 
quire  <?e(0)  =  <7ot(0)  =  gvd(0)  =  0  so  that  a  zero  strain  will  yield  a  zero  stress  according  to  our 
stress-strain  relation.  Based  on  additional  studies  detailed  in  [22],  we  chose  an  exponential 
form  Y(s)  —  e~r  for  the  memory  kernel.  Such  an  exponential  form  generates  totally  nested 
hysteresis  loops  in  the  stress-strain  curves,  a  feature  also  observed  in  our  data. 

The  parameters  r,  Ei:  a*  and  6,  were  estimated  by  setting  up  a  least  squares  minimization 
problem  using  one  or  two  of  the  outer  stress  strain  data  loops  from  the  experiments.  Figure 
3  shows  results  for  a  highly  filled  sample  in  Type  I  and  Type  II  experiments,  respectively. 
While  the  inverse  problem  was  “trained”  on  the  largest  loop  only  in  Type  I,  and  the  two  outer 
loops  in  the  Type  II  experiment,  we  can  see  that  there  is  a  very  good  agreement  between  the 
data  and  the  model  predicted  inside  loops.  Similar  results  were  obtained  for  medium-filled 
natural  rubber  and  silica- filled  silicon  samples.  More  details  are  given  in  [21]. 

Our  hysteretic  model  (11)-(12)  was  also  tested  on  a  series  of  dynamical  experiments  with 
different  types  of  filled  rubber  samples.  We  used  the  same  free  release  experiments  that 
were  described  earlier.  Since  the  dynamical  behavior  of  the  unfilled  natural  rubber  sample 
was  adequately  described  by  the  basic  model  without  hysteresis,  we  began  our  hysteresis 
investigations  using  the  lightly  filled  sample.  For  given  p,  CD:  ge ,  gv  and  Y  we  solve  the  partial 
differential  equation  (11)-(12)  forward  in  time,  and  obtain  the  displacement  u(t,  x),  0  <  x  < 
£.  The  data  collected  in  these  experiments  provides  us  with  the  force  at  the  top  of  the  sample 
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Figure  3:  (a)  Model  prediction  for  Type  I  CB-h  data  and  (b)  Model  prediction  for  Type  II 
CB-h  data. 


( x  =  0),  and  we  compare  it  to  the  force  predicted  by  the  model  at  the  same  point,  given  by 

AMt,  0)  =  A  0)  +  M^(t,  <•))  + 1‘  Y(t  -  0).  |£(»,  "))*  . 

where  we  use  our  computed  solution  u(t,x)  to  find  |^(f,  0)  and  j^(f,  0).  Our  goal  is 
to  find  p,CD,ge,gv  and  Y  so  that  the  model  predicted  force  at  x  —  0  best  matches  the 
data  collected  by  the  load  cell.  A  parameter  identification  problem  was  set  up  to  find 
p,  CD,  Ei,  E2,  E3,  ai,  a2,  a3,  b±,  b2,  63,  r  (collectively  denoted  by  q)  such  that 

N 

J(q)  =  J2  I zi  -  Acv(ti,0-,q)\2 

i= 1 

is  minimized.  Here,  once  again,  the  zi:  i  =  1  IV,  are  the  data  collected  by  the  load 
cell,  and  <r(fj,0;g)  is  given  as  described  above.  The  particular  forms  for  ge,gv  and  Y  were 
motivated  by  their  success  in  the  quasi-static  case.  In  our  computations  we  used  linear  splines 
for  spatial  discretization.  In  solving  the  system  (11)-(12)  forward  in  time,  the  treatment  of 
the  hysteresis  integral  proved  to  be  very  time  consuming.  Since  this  computation  needs  to 
be  repeated  for  each  set  of  parameters  during  the  optimization,  the  time  required  for  the 
computational  parameter  identification  process  was  very  extensive.  Hence,  we  formulated 
an  equivalent  system  to  (11)-(12)  using  an  internal  variable  £1  =  e\ (i;r),  and  used  it  in  the 
above  framework  for  our  subsequent  calculations.  This  system  is  given  by 

inV'  (13) 

1  d  du  d2u 

£l  =  _7  ei  +  Tt9’(di’Mdi) 


(14) 
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Figure  4:  Fit  to  data  with  model  using  a  nonlinear  ge,  gv  :  (a)  Time  domain  approximation 
and  (b)  the  power  spectrum  of  the  approximation  and  the  data  for  a  lightly  filled  sample 
with  3  lb  tip  mass. 


u(t,0)=0,  u(0,x)  =  ipo,  —  (0,z)  =  0,  £i(0,x)  =  0,  (15) 

where,  in  general,  the  derivatives  of  gv  in  (14)  are  distributional  in  the  sense  described  earlier. 
The  parameter  identification  problem  was  solved  using  MATLAB  optimization  routines. 
Typical  results  for  the  identification  problem  are  shown  in  Figure  4,  where  a  3  lb  tip  mass 
was  used  at  the  bottom  of  the  sample.  The  identification  problem  was  also  run  on  data 
obtained  from  experiments  with  2  lb  and  1  lb  tip  mass  added  to  the  lightly  filled  sample. 
We  found  a  similar  good  agreement  between  the  data  and  the  model  [21],  However,  the 
identified  coefficients  were  not  entirely  consistent  across  experiments,  although  they  should 
describe  the  same  material.  This  variation  is  probably  caused  by  the  substantially  different 
stress  and  strain  rate  ranges  involved  in  the  experiments,  and  suggests  that  the  model  should 
be  refined  to  account  for  these  differences. 

We  repeated  the  experiment  with  a  highly  filled  sample  having  a  9.29  lb  tip  mass.  In  this 
case  our  best  fit  depicted  in  Figure  5  (a)  has  deteriorated.  Thus  we  turned  to  a  modification 
of  our  model  to  obtain  a  better  approximation  to  the  force  data. 

This  modification  was  motivated  by  our  earlier  efforts  as  well  as  those  reported  in  the 
literature  on  similar  problems  involving  hysteretic  effects  and  internal  strain  models  (for 
details,  see  the  discussions  in  [21]).  We  next  considered  the  model  using  two  internal  variables 
Ti),  £2(t]  t2)  with  different  decay  parameters  ti,t2,  given  by 


„d2u  d3u 

PW~  cCDem? 

1  d  . 

£i  — - £i  +  ~r:9v\ 

Ti  dt 


d_ 

dx 


*  ,du  N 

Acge{—)  +  Ace  i  +  Ac£2 


du  d2u  . 
dx1  dtdx 


=  m 


in  V 


(16) 


(17) 
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highly  filled  sample  highly  filled  rubber  sample 


Figure  5:  (a)  Best  approximation  for  the  highly  filled  sample  with  one  internal  variable  using 
data  in  the  time  interval  [0,  0.95]  and  (b)  Best  approximation  for  the  highly  filled  sample 
with  two  internal  variables  on  the  full  time  interval.  Identification  based  on  data  from  the 
time  interval  [0,0.95]. 


1  d  ,du  d2u  , 
2  r2  2  +  dtPv  dx:  dtdx' 


(18) 


u(t,0)  =  0,  u(0,x)  =  <p0,  —  (0,z)  =  0,  £i(0,  x)  =  0,  e2(0,x)=0.  (19) 

After  parameter  estimation,  this  model  provided  a  very  satisfactory  fit  to  the  data  with 
the  relative  error  between  the  data  and  the  computed  force  being  4.1%.  This  fit,  depicted  in 
Figure  5  (b),  was  obtained  by  using  data  only  from  the  time  interval  [0,  0.95],  and  resulted 
in  a  model  that  accurately  simulated  the  data  on  the  interval  [0, 1.6]. 


2.3  Nonlinear  internal  variable  models 


Our  last  result  for  the  highly  filled  sample  and  the  variation  in  the  parameters  for  the  lightly 
filled  rubber  rod  suggested  that  we  might  try  to  generalize  our  model  to  better  describe 
the  behavior  of  highly  hysteretic  samples.  Thus  we  next  considered  internal  variable  models 
with  nonlinear  internal  dynamics 


„d2u  A  ri  d3u 

PW~  cCDdtd^~ 

£i  =  ~9e{e  l)  +  jt9v{ 


dx  [AMdx)  +  Ac£\ 
hi  d2u  , 


dx ’  dtdx ' 


=  F(t ) 


in  V* 


u(t,  0)  =  0,  u(0,  x)  =  ipo,  —  (0,x)  =  0,  £i(0,x)=0. 

In  [23]  we  showed  that  this  system  has  a  unique  weak  solution.  We  also  note  that  we  in  fact 
generalized  the  previous  existence-uniqueness  result  in  the  sense  that  we  no  longer  require 
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monotonicity  of  the  nonlinear  functions  ge-,  9vu  9vd-  Instead,  they  are  assumed  to  satisfy  a 
local  Lipschitz  property  and  other  assumptions  similar  to  those  used  before.  We  remark 
that  similar  techniques  were  successfully  employed  to  establish  existence-uniqueness  of  weak 
solutions  for  linear  evolution  equations  of  second  order  in  t  in  [24]  and  for  semilinear  second 
order  evolution  equations  where  the  nonlinear  term  satisfies  a  global  Lipschitz  condition 
in  [25].  In  [26]  we  used  similar  techniques  to  study  a  nonlinear  beam  equation  where  the 
nonlinearity  satisfies  only  a  local  Lipschitz  condition. 


2.4  Model  development,  experimental  and  numerical  results  for  simple  shear 

The  previous  model  development  and  identification  process  was  also  carried  out  for  filled 
elastomers  undergoing  simple  shear  deformation.  A  dynamic  experiment  was  designed  that 
involved  a  “double  sandwich”  fixture  with  layers  of  filled  rubber  at  the  interfaces  (see  Figure 
6).  The  side  bars  were  fixed,  while  the  middle  bar  was  perturbed  by  an  impulsive  hammer 
hit,  and  accelerometer  and  load  cell  data  were  collected  as  in  the  experiments  for  extension. 
For  simple  shear  the  model  developed  and  used  was 


pA, 


d2u  d 


dt2  dx 
d2u 


a  d2u 

+ AfiDm~x + A£l 


,du  d2u 

9’(di)  +  CDdW~x+ei 


+  fit)  +  Mg 


d  1  d  ,du  d2u  . 

df£l  r£l  dt^v  dx:  dtdx 

du 

u(t,  0)  =  0,  u(0,x)  =  (po,  —(0,x)  =  0,  £i(0,  x)  =  0. 


Load  cell 


extra  mass 


hammer  hit 


Figure  6:  Schematic  of  the  initial  shear  experimental  device. 
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In  the  identification  process  we  used  the  accelerometer  data  to  approximate  the  unknown 
parameters,  i.e. ,  we  solved  the  least-square  minimization  problem 

1  N  B2v 

Aq)  =  g  -  q 

where  a,  represent  the  collected  data,  and  q  is  a  vector  of  unknown  coefficients.  The  ex¬ 
periment  was  performed  with  lightly  filled  and  highly  filled  rubber  samples  with  different 
amounts  of  extra  mass.  In  general,  we  found  that  the  data  collected  with  no  extra  mass 
on  the  sample  could  be  approximated  well  with  the  above  model  even  when  the  viscoelastic 
response  gv  is  assumed  to  be  linear.  This  is  not  surprising,  since  the  maximum  strain  levels 
in  these  type  of  experiments  were  below  10%.  Our  best  fit  for  the  highly  filled  sample  with 
no  extra  mass  is  depicted  in  Figure  7.  However,  when  extra  mass  was  involved  to  achieve 
larger  deformations  and  strain  levels,  the  experiment  did  not  provide  suitable  data,  since 
additional  “tilting”  modes  became  excited  that  could  not  be  accounted  for  by  the  above 
one-dimensional  model. 


A225  (b)  1 


time 


Figure  7:  Highly  filled  rubber  sample  in  shear  with  no  extra  mass 


The  experiment  was  redesigned  to  enforce  simple  shear  in  the  sample  with  no  additional 
deformations.  These  new  experiments  involved  four  “double  sandwich”  samples  with  highly 
filled  rubber  at  the  end  of  four  arms  of  a  fixture.  The  arms  were  latched  down  to  achieve 
a  prescribed  initial  strain  (50%-100%)  in  the  sample  and  subsequently  released.  We  again 
collected  accelerometer  and  load  cell  data.  The  identification  problem  was  first  performed 
for  the  data  obtained  for  100%  initial  strain.  We  tried  different  hysteretic  stress-strain  rela¬ 
tionships  in  our  numerical  simulations.  We  note  that  a  linear  viscoelastic  response  function 
gv  was  no  longer  adequate  to  describe  the  data  in  this  higher  strain  regime.  Thus,  we  initially 
assumed 


,du.  d2u  .  1  ,du. 

a  =  9€{di)  +  CDdiai  +  £'’  Sl  +  7Sl=9’(&)’ 
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where  ge  and  gv  are  cubic  polynomials  and  gv  does  not  depend  on  i,  (i.e. ,  gvi  =  gv(i  in  (9)). 
The  parameters  that  provided  the  best  fit  in  this  case  did  not  work  well  for  the  70%  and 
50%  initial  strain  experiments.  Thus  we  next  used  the  full  nonlinear  hysteretic  model  (i.e., 
9vi  %  9vd  in  (9))  to  approximate  the  100%  initial  strain  data.  The  best  fit  in  this  case 
is  depicted  in  Figure  8  (a).  We  found  that  the  set  of  parameters  identified  for  this  case 
described  the  data  with  70%  initial  strain  with  very  good  accuracy  (see  Figure  8  (b)). 


Acceleration  data  and  the  approximated  values 


Acceleration  data  and  the  approximated  values 


Figure  8:  (a)  Approximation  result  for  the  highly  filled  sample  in  shear  with  ~  100%  strain 
initial  condition  and  (b)  The  same  parameter  set  simulating  the  highly  filled  sample  in  shear 
with  70  %  strain  initial  condition. 


2.5  Molecular  based  reptation  models 

The  above  models  and  the  nonlinear  extensions  of  Section  2.3  have  not,  to  date,  provided 
insight  into  the  underlying  mechanisms  for  tensile  and/or  shear  deformations  in  filled  rub¬ 
ber.  This  is  not  unexpected  since  the  approaches  described  above  are  based  on  pseudo- 
phenomenological  formulations.  More  recently  ([5,  27]),  we  have  turned  to  a  different  ap¬ 
proach  based  on  molecular  arguments  which,  as  we  shall  see,  lead  precisely  to  the  class  of 
models  based  on  a  Boltzmann  hysteresis  formulation.  As  usual,  one  begins  with  force  and 
moment  balance  and  seeks  constitutive  laws  (such  as  (10))  for  the  viscoelastic  stress  term 

O vis co  m 

where  £  =  is  the  infinitesimal  strain  and  si  is  an  “internal  strain”  variable  on  which  avisco 
depends  in  an  hysteretic  manner.  As  described  above,  we  found  that  a  reasonable  description 
of  the  data  could  be  given  with  the  typical  stress-strain  relationship 

<r(*)  =  9e(e(t),e(t))  +  Jo  7  e-L^—gv(e(s),i(s))ds, 
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where  r  is  a  relaxation  parameter,  gv  is  as  defined  in  (9),  with  cubic  polynomials  gvi,gvd,  and 
9e  =  9e( §y)  +  Cd-§^.  as  in  (4).  We  have  already  observed  that  this  expression  is  equivalent 
to 


a(t)  =  ge(e(t),  s(t))  +  ysrft;  r), 

where,  for  a  given  “relaxation  parameter”  r,  the  internal  strain  £1(t;r)  satisfies  (8).  In  fact, 
we  found  that  highly  filled  rubbers  required  multiple  relaxation  times  Ti,t2  as  in  (16)-(19) 
to  obtain  good  model  fits  to  the  data.  As  we  shall  describe,  molecular  based  formulations, 
where  microscopic  relaxation  parameters  vary  across  the  population  of  molecules  in  the 
material,  lead  to  internal  dynamics  of  the  form  (8), (9).  When  combined  with  a  Prohorov 
metric  framework  (see  [28,  29])  for  uncertainty  in  internal  dynamics,  these  ideas  lead  to  the 
computational  models  we  have  used  above.  Indeed,  the  molecular  based  approach  leads  to 
a  general  class  of  models  with  uncertainty  or  randomness  in  the  stress 

a(t,x;P)  =  ge(e(t,x),i(t,x))  +  7  jT_ei(i,a;;r)dP(r),  (20) 

where  P  is  a  probability  distribution  over  the  set  T  of  possible  relaxation  parameters,  and 
Si(t;  t)  satisfies,  for  each  r  G  T, 

£i(f,  x]  r)  H — £1  (t,  x ;  r)  =  i(t,  x)h(e(t ,  x)). 

T 

For  the  reptation  model  derivation  in  [5],  one  begins  with  the  Doi/Edwards  [3]  stick-slip 
molecular  models  as  embodied  in  the  continuous  tube  reptation  models  of  Johnson/Stacer 
[4]  wherein  polymer  materials  such  as  rubber  are  postulated  to  be  composed  of  two  types  of 
molecules.  In  tensile  deformations,  one  denotes  by  L(t)  the  length  of  chemically  cross-linked 
or  CC  molecules,  while  l(t)  denotes  the  length  of  physically  constrained  or  PC  molecules. 
To  use  stick-slip  models  in  continuum  simulations  of  reptation  in  rubbers,  one  considers 
networks  of  “cells”  or  boxes  of  parallel-sided  CC  boxes  and  PC  boxes  with  sides  of  length 
(principal  stretches) 

>  duc  ,  „  _  duv 

Ac  —  l  +  £  —  1  +  — — ,  Xp  —  1  +  £\  —  1  +  — — , 
ox  ox 

respectively.  Here  uc  denotes  the  deformations  of  the  CC  box  and  up  denotes  the  deforma¬ 
tions  of  the  PC  box.  Using  a  linear  stick-slip  assumption  as  in  [4],  and  strain  energy  densities 
based  on  experiments  of  Young  and  Danik  (see  [5,  27]  for  details),  one  obtains  as  a  limit  of 
PC  response  to  step  tensile  deformations  of  the  CC  molecules,  the  s,  £1  coupled  dynamics 

1  _  .  1  +  £1 

£l  H - £\  —  £— — ■ — . 

r  1  +  £ 

However,  if  one  replaces  the  linear  assumption  of  [27]  by  a  nonlinear  stick-slip  hypothesis 
(which  is  the  basis  of  the  work  in  [5]),  one  obtains  a  more  general  nonlinear,  dynamical 
relationship  between  £  and  £1  given  by 

£i-| — £1  =  £/((  1  -F  £i)/(l  H-  £))- 
r 
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Expansion  and  truncation  of  higher  order  terms  lead  to  equations  of  the  form 

£ l  H - £\  —  £'(o!o  +  OL\£  +  Oi2^2  +  Q!3£3),  (21) 

r 

which  are  of  the  same  form  as  the  internal  variable  model  (8), (9)  with  gvi  and  gvd  cubic 
polynomials.  For  the  corresponding  contributions  to  a  from  the  strain  energy  densities 
of  Young-Danik/  Johnson-Stacer  with  the  nonlinear  stick-slip  hypothesis,  one  obtains  a 
contribution  to  the  rate  independent  strain  gsv  (after  expanding  /  in  a  Taylor  series  and 
dropping  higher  order  terms)  of  the  form 

gl(e ,  £i)  =  9cubic(£ )  +  7i£ij 

where  £i  is  as  before  (i.e. ,  the  internal  strain  satisfying  (21)).  Thus,  the  total  stress-strain 
relationship  can  be  written  in  the  form  (20).  If  the  measure  P  of  (20)  has  atoms  at  T\  and 
r2,  (i.e.,  the  measure  is  composed  of  Dirac  measures  concentrated  at  iq  and  r2),  then  the 
constitutive  law  leads  precisely  to  the  model 

cr(t,  X]  P )  =  ge(e(t ,  x ),  i(t,  x))  +  7i£i(t,  x;  rx)  +  7 2£2(t,  x;  r2), 

which  was  used  in  the  data  fit  of  Section  2.2  with  model  (16)-(19). 

3  Concluding  Remarks 

In  the  above  note  we  outlined  our  progress  to  date  in  the  development  of  nonlinear  dynamic 
models  for  inactive  filled  elastomers.  Substantial  experimental  validation  for  our  approach  is 
provided  both  in  the  quasi-static  and  dynamic  cases  in  uniaxial  tension  and  in  the  dynamic 
case  in  simple  shear.  Current  efforts  involve  refinements  to  these  models  and  a  compar¬ 
ison  with  newly  developed  molecular  based  models  [5,  27]  that  emphasize  uncertainty  or 
randomness  across  populations  of  molecules  in  a  heterogeneous  material. 
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